Diversity enabling equilibration: disorder and the ground state in artificial spin ice 
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We report a novel approach to the question of whether and how the ground state can be achieved 
in square artificial spin ices where frustration is incomplete. We identify two sources of randomness 
that affect the approach to ground state: quenched disorder in the island response to fields and 
randomness in the sequence of driving fields. Numerical simulations show that quenched disorder 
can lead to final states with lower energy, and randomness in the sequence of driving fields always 
lowers the final energy attained by the system. We use a network picture to understand these two 
effects: disorder in island responses creates new dynamical pathways, and a random sequence of 
driving fields allows more pathways to be followed. 

PACS numbers: 75.50.Lk, 75.75.-c, 75.78.-n 



The ability to control artificial spin ice [l|-[15j and re- 
lated systems [H, [I?} has made them a useful testing 
ground for many types of physics. For example, the 
effective temperature formalism used for granular ma- 
terials can be studied in the context of artificial spin 
ice [1, [TH, EH- The submicron Ising-like magnetic is- 
lands of artificial spin ice are arranged so the dipolar 
inter-island interactions are frustrated, leading to a com- 
plex energy landscape with many states nearly degener- 
ate. In square artificial spin ices, the interactions within 
each vertex of the lattice of islands are inequivalent, lead- 
ing to a well-defined, two-fold degenerate ground state, 
as shown in Fig. []Ja). Large domains of ground state 
ordering have been observed in samples that have under- 
gone thermal annealing during growth Q. However, once 
grown, the islands are large enough to be athermal and 
dynamics are driven entirely by external fields. 

One difference between athermal driven and thermal 
dynamics is that an external driving field acts uniformly 
on all islands, whereas thermally driven individual is- 
lands may reverse stochastically. A consequence of hav- 
ing global and deterministic driving - rather than local- 
ized and random - is that dynamics are constrained and 
even though the inequivalence of interactions lifts degen- 
eracies in the energy landscape, the field driven dynamics 
cannot neccessarily evolve the ice to a low energy state. 
Experimentally, a sequence of fields applied to a satu- 
rated configuration (shown in Fig.QJb)) gives a final state 
with only short range ground state ordering [l], 0, H, [lo[ . 

Although the driving field is global, a spin's behavior 
is also determined by local factors: its interactions with 
other spins and its (intrinsic) switching field. Disorder in 
these affects the response to driving. Disorder has previ- 
ously been shown experimentally to be important for the 
response of square artificial spin ice to fields [9] but until 
now, the mechanisms have been unclear. Here, we focus 
on the simplest type of disorder, which is to allow the 
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FIG. 1. The array geometry studied, shown for 6 x 6 is- 
lands, (a) The ground state configuration in which spins are 
arranged in a microvortex configuration. The ground state is 
two- fold degenerate, with a global spin flip giving the other 
ground state configuration, (b) A saturated configuration, in 
which all spins have a positive projection onto an axis. There 
are four saturated configurations. 



switching fields to vary in magnitude from island to is- 
land, as would result from roughness in the islands. This 
variation in switching fields does not modify the inter- 
actions, and frustration and degeneracies are unaffected. 
Instead, it modifies the way the system explores its en- 
ergy landscape. We study the effects of disorder using 
numerical simulations which give full control over island 
properties, and by mapping transitions onto a directed 
network which can be studied using network theoretic 
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[22|. This latter approach provides a concep- 
tual framework which can be applied to other systems. 

We first describe the effect of switching field disorder 
on the final energy of the system, when a uniformly rotat- 
ing protocol is used. Then, we interpret the results with 
the support of the network description, which is also the 
natural framework to study changes in the field protocol. 
Finally, we study field protocols with a random sequence 
of field angles, and we show they allow the system to 
achieve lower-energy final states than can occur under 
rotating protocols, in a way that is more robust against 
quenched disorder. 

In our simulations, the magnetic islands are treated 
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as dimensionless Ising spins Si (S = 1) that interact as 
point dipoles. The energy of spin i is given by its Zeeman 
energy, —h-Si, and the sum of dipolar interactions with 
all other spins: 
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spin z can flip if 
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where h[^ t = h + fe^ and ft^ is the island's intrinsic 



switching field. 

Evolution under an external field h proceeds by ran- 
domly selecting a spin satisfying criterion ([2]), flipping it, 
then checking the switching condition again for all spins, 
until no further flips are possible. We denote the tran- 
sition under a field h from initial state i to final state / 
by i —> f • (More than one final state may be possible; 
alternatively / may be the same state as i if no spins can 
flip.) The transition i — >• / is the building block of the 
simulated dynamics and the network picture, discussed 
below. 

In the absence of disorder, all islands have the same 
switching field. We implement disorder by drawing the 
hc^ from a Gaussian distribution characterized by stan- 
dard deviation a. We work in units where l/(47r/xo) an d 
the nearest-neighbor distance are set to unity, and the 
mean h c value is 11.25. 

In a perfect system with edge geometry as shown in 
Fig. [TJ a rotating field with constant amplitude can in- 
duce non-trivial dynamics only in a narrow range of field 
amplitudes, Ah ~ 2. This is shown for a 20 x 20 spin 
system as blue circles in Fig. O As discussed in Ref. [6|, 
smaller fields are unable to flip spins and larger fields 
force the magnetization to track the field. Dynamics al- 
ways start at array edges where spins have fewer neigh- 
bors. The minimum of Edip(h) at h ~ 10.5 is due to a 
very regular and specific process of spin flips that "in- 
vade" from the array edges. 

Let us now consider the effect of switching field disor- 
der, with <j = 1.875. This disorder strength is consistent 
with exp erimental data for the square lattices studied 

Dis- 
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here [23] and disconnected kagome ices 
order allows dynamics to start inside the array at sites 
where "loose" spins with smaller switching fields are lo- 
cated. Similar disorder effects have been seen in simula- 
tions of vortices in nanopatterned superconductors jTjj. 
The uniformly rotating field protocol now produces the 
curve Edi p (h) given by inverted triangles in Fig. El The 
increase in Ah indicates that the available configuration 
space is enlarged: dynamical processes forbidden in the 
perfect system - such as the flipping of isolated spins in 
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FIG. 2. Disorder in island switching fields or the sequence 
of applied fields allows the system to reach final states with 
low energy. As indicated by the legend, symbols represent 
rotating and random field protocols acting on perfect and 
disordered systems. Averages are made over 20 disorder real- 
izations. Error bars represent one standard deviation. 



the array bulk - are allowed when disorder is present. 
However, disorder can also destroy transitions. For ex- 
ample, the orderly process of spin flips at h ~ 10.5 in the 
perfect system cannot occur in the disordered system. 
This leads to higher energy at these fields. 

Up to now, the discussion has focussed on the dynamics 
of spins in the array. However, an alternative viewpoint 
- alluded to in the discussion of transitions i —> f - is 
to see the action of the applied field as "transporting" 
the system from one spin configuration state to another, 
in a path through the space of all configurations. This 
picture can be interpreted as a network, in which net- 
work nodes are spin configurations and a directed link 
(i —> f) exists if applying some field to configuration i 
allows it to transition to configuration /, according to 
the dynamics described above. The advantage of this 
approach is that it allows us to discuss dynamics in a 
way that is not restricted to a particular field protocol. 
Related approaches to mapping dynamics onto networks 
have proved fruitful in the study of other geometrically 
frustrated systems [jj 



25l | and the random field Ising 
model [26|-l28|. 

We restrict our considerations to fields of strength 
h = 11.5. Using a different field amplitude would give 
a different network. In principle, the field can take any 
angle in (0,27r] but we consider only the discrete set of 
field angles 6 = n7r/128, n = 1,2,... 256. As seen in the 
Supplementary Material [29|, this set of field angles can 
be expected to give results close to the limit of continuous 
0. For each configuration (node), transitions are calcu- 
lated at all field angles 0. A field protocol {0i,02>---} 
corresponds to a path on the network where at each step 
of the path only a link corresponding to the angle 0i can 
be followed. 

Since the number of configuration states increases ex- 
ponentially with the size of the lattice, we consider a 
small, feasible to analyze, 4x4 array, with 16 islands 
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FIG. 3. The subnetwork of states that can be accessed from 
the +x saturated configuration for (a) a perfect system (where 
the saturated configuration is the large red node) and (b) a 
typical disordered system. 



and 2 16 total configurations. We are able to obtain the 
exact network representation of the dynamics of such a 
system, via enumeration. Although the dynamics of the 
4x4 system are not exactly the same as those of the 
20 x 20 system, they are governed by the same basic rules, 
and results for the smaller system are applicable to the 
larger one [29|. The networks are stored as adjacency 
matrices. Aif = 1 if the transition i —> f is allowed for 
some Aif = otherwise. These are sparse matrices, 
with ~ 10 6 non-zero entries in a 2 16 x 2 16 matrix. Disor- 
der is implemented in a similar way to the simulations. 
However, because of the small number of islands, the ac- 
tual mean of the h c values drawn from the distribution 
is 11.25 + A, where A can be large enough to affect the 
network properties. We subtract A from each h c value to 
fix the mean at 11.25. Each disorder realization gives a 
different network, but we find similar network properties 
for disorder realizations with similar standard deviations 
in h c . 

Fig. [3ja) shows the subnetwork of states reachable 
from the +x saturated configuration in a perfect system. 
All the links that exist for h = 11.5, = n7r/128 are 
shown. The most striking feature of the subnetwork is 
that it is very limited, containing only five nodes. The 
saturated configuration (indicated by the large red node) 
can evolve into one of two states, each of which can evolve 
into a single state, giving two possible final states for 
the system. Self-loops indicate that certain field angles 
are unable to induce transitions. There are no "return" 
paths - transitions are irreversible. Fig. [3fb) shows the 
subnetwork of states reachable from the same node, for 
a typical disorder realization. The subnetwork contains 
1814 nodes, three orders of magnitude more than the 
perfect system. In this subnetwork, there are paths into 
the saturated configuration, as well as out of it. Other 
realizations of disorder give similar results, with the num- 
ber of configurations reachable from the saturated state 
ranging from 1758 to 12516 in a study of ten disorder 
realizations. 

Turning to the network as a whole, we see that disorder 
"re-wires" the network significantly, increasing the total 
number of links from 736,720 to 949,216. This is not 
simply the addition of links: disorder removes 358,934 



links from the perfect system, and creates 571,430 links. 
As already mentioned, one result of this re-wiring is an 
increase in the reversibility of dynamical pathways. This 
can be described quantitatively using strongly connected 
components (SCCs). In a directed network, if nodes A 
and B are in the same SCC, there exists a path from A 
to B and vice-versa. The SCCs are determined using the 
algorithm in Ref. [30j |. as implemented by the software 
package Mathematica. Fig. HJa) shows the distributions 
of SCC sizes for perfect and disordered systems. In the 
perfect system, the largest SCC contains 3 nodes and 
the saturated configurations are in SCCs of size 1. In the 
disordered system the largest SCC contains 1628 nodes, 
and includes the saturated configurations. In terms of dy- 
namics, this means that in the disordered system there 
exists a sequence of fields to take a saturated configu- 
ration to any one of 1627 configurations, and also a se- 
quence back to saturation. Other disorder realizations 
give maximum SCC sizes ranging from 606 to 12386. 

On the other hand, as seen in Fig.HJb), the in- and out- 
degree distributions are not affected much by disorder. In 
both perfect and disordered systems, most nodes have a 
low in-degree, indicating that most configuration states 
can only be reached from a small number of other states. 
The peak in the out-degree distribution is at fc out ~ 10, 
with only a few nodes having very low out-degree (stable 
configurations) or very high out degree (unstable config- 
urations). As seen in the insets to Fig. HJb), nodes that 
have low (high) degree in the perfect system typically also 
have low (high) degree in the disordered system. The de- 
gree is increased more often than it is decreased because 
disorder adds links to the system which are distributed 
among the nodes. 

While diversity in switching fields can open new links 
between states, a link i — >• / is active only for certain 
orientations of the field, Of^j < < Of^j. If the system 
is in state z, but the applied field is not at an angle in the 
active interval, the link is not "operative". In this way, 
the field protocol is essential in determining dynamics. 
Changing the protocol can dramatically modify the final 
configurations attained. 

We now show evidence of the importance of field pro- 
tocol by simulating alternative protocols. The first ex- 
ample, examined in Fig. [5j takes advantage of the way 
disorder breaks symmetry between clockwise (CW) and 
anticlockwise (ACW) rotations. In the perfect system, a 
uniform rotation protocol gives the same result whether 
the rotation is CW or ACW, because of reflection sym- 
metry. When disorder is introduced, the two senses of 
rotation are no longer equivalent for any particular re- 
alization of disorder, but symmetry is restored when we 
average over disorder, as seen in Fig. [5] We can go a step 
further and consider a protocol that samples both senses 
of rotation by repeatedly making a number of rotations 
CW and then ACW. In the perfect system, this alter- 
nating CW/ACW protocol is equivalent to one with a 
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FIG. 4. (a) The distributions of strongly connected compo- 
nent sizes for the perfect (orange/light gray squares) and dis- 
ordered (blue/dark gray circles) systems, (b) The in- and 
out-degree distributions for the perfect (orange/light gray 
squares) and disordered (blue/dark gray circles) systems. In- 
sets: The in- and out-degree of 1000 randomly chosen configu- 
rations, in the disordered vs the perfect system. The red line 
separates those nodes whose degree is increased from those 
whose degree is decreased by disorder. 



single sense of rotation. In a disordered system, it leads 
to lower energy states with greater ground-state order- 
ing. The key point is that this simple modification of the 
uniform protocol, and the presence of disorder, allow the 
system to avoid trapping by higher-energy local minima. 

The alternating CW/ACW protocol breaks the regular 
sequence of applied fields to allow the system to traverse 
regions of its phase space network inaccessible to the uni- 
form rotating protocol. Even more effective is a random 
sequence of field orientations, corresponding to a random 
walk on the network. The randomness introduced by this 
protocol is extrinsic and controllable, unlike the random- 
ness introduced by switching field disorder. Results for 
random protocols are shown in Fig. [2j The diamonds 
show results for a perfect ice, and the squares show re- 
sults for a disordered ice. In both cases the final energy 
is lower than that reached by a uniformly rotating pro- 
tocol in the same system, and the fraction of vertices in 
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FIG. 5. For an undisordered system, clockwise rotating fields, 
anticlockwise rotating fields, and a field protocol that alter- 
nates sense of rotation every 2 cycles all give the same final 
energy. When disorder is included, the clockwise and anti- 
clockwise rotating fields give the same results on average, but 
the alternating field protocol leads to a lower final energy. 
Symbols for each protocol and system are indicated in the 
legend. Averages are made over 100 disorder realizations and 
error bars represent one standard deviation. 



a ground-state configuration is higher (not shown). The 
large error bars for h > 11.5 indicate the spread in fi- 
nal configurations attained, corresponding to a variety 
of pathways taken on the underlying dynamics network. 
Indeed, for h > 13.5, the simulations do not attain a sin- 
gle final state, and the averages of Fig. [2] are made over 
configurations attained after the arbitrary time of 2000 
field applications. 

It is interesting to contrast these results for an ather- 
mal ice to those presented recently for thermal ices. In 
Ref. [7[ , experimentally observed large scale ordering was 
attributed to thermally driven processes occurring dur- 
ing growth stages where the islands were small enough to 
undergo thermal fluctuations. Thermal fluctuations are 
often thought of as stochastic local fields. In the athermal 
system considered here, they are replaced by either a ran- 
dom sequence of global fields or frozen randomness in the 
form of a distribution of switching fields. The unifying 
principle for these two types of randomness is the con- 
cept of links between configurations. The links followed 
depend on both the switching fields and the driving field. 
Diversity in either of these enables a wider exploration 
of configuration space. 

This ability to achieve long range ground state order- 
ing in the athermal system suggests a process similar in 
certain respects to the "order by disorder" proposed by 
Villain for equilibrium systems [31]. In order by disorder 
phenomena, disorder lifts degeneracies in frustrated sys- 
tems and a ground state, inaccessible in a pure system at 
zero temperature can be acheieved by introducing either 
quenched or thermal disorder. This is an equilibrium 
phenomenon where fluctuations select a specific ground 
state from an ensemble of degenerate states. By way 
of contrast, our athermal nonequilibrium system with- 
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out disorder has a single ground state and strongly con- 
strained dynamics. While disorder does not make the 
ground state accessible [23|, it does remove constraints 
by adding links between states. With the correct field 
protocol, these links can be followed, allowing the sys- 
tem to relax into a state with a high level of ground state 
ordering. 
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